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We study the kinetic theory of driven granular gases, taking into account both translational and 
rotational degrees of freedom. We obtain the high-energy tail of the stationary bivariate energy 
distribution, depending on the total energy E and the ratio x — ^/E w / E of rotational energy E w to 
total energy. Extremely energetic particles have a unique and well-defined distribution f(x) which 
has several remarkable features: x is not uniformly distributed as in molecular gases; f(x) is not 
smooth but has multiple singularities. The latter behavior is sensitive to material properties such 
as the collision parameters, the moment of inertia and the collision rate. Interestingly, there are 
preferred ratios of rotational-to-total energy. In general, f(x) is strongly correlated with energy 
. and the deviations from a uniform distribution grow with energy. We also solve for the energy 

Oh' distribution of freely cooling Maxwell Molecules and find qualitatively similar behavior. 
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I. INTRODUCTION 

Energy dissipation has profound consequences in granular materials, especially in dilute gases, where the dynamics 
are controlled by collisions [H, 0, H[ . Dissipation is responsible for many interesting collective phenomena including 
clustering [J, 0, 0, 0, 0] , formation of shocks d, E3, El, EH EH , and hydrodynamic instabilities [HI, ElL A nother 
co nseq uence is the anomalous statistical physics that includes the non-Maxwellian velocity distributions 1(| E3, El, 
El, US El and the breakdown of energy equipartition in mixtures |22l . [23[ . 

For an elastic gas in equilibrium, the temperature, defined as the average kinetic energy, characterizes the entire 
distribution function including all of the moments, the bulk of the distribution, as well as the tail of the distribution. 
Outside of equilibrium, the temperature is not sufficient to characterize the energy distribution. Granular gases are 
I , inherently out of equilibrium and a complete characterization must therefore include the behavior of typical particles, 

■ the behavior of energetic particles, as well as the moments of the distribution. For example, the energy distribution 
may have power-law tails with divergent high-order moments [2~3 . [25l [26| and consequently, the moments exhibit 

■ multiscaling [27| . Generally, nonequilibrium effects are pronounced in the absence of energy input to balance the 
dissipation but can be suppressed by injection of energy where the deviation from a Maxwellian distribution affects 

. only extremely energetic particles [l?], HI, HtJ H(| • 

While there is substantial understan ding of the energy distribution of frictionless granular gases, much less is known 
theoretically [H], H2, HI, H3, HI, HI, H3, l38| | and experimentally [39L |40L l4lj when the rotational degrees of freedom are 
taken into account. It is difficult to measure the rotational motion experimentally, and the few available measurements 
are restricted to two-dimensions. Surface roughness and friction have important consequences and the hydrodynamic 
theory [H, HI, H3, |4f| must be modified, if the particles have spin [4(| . Equipartition does not hold for the average 
rotational and translational temperature neither in the free cooling case [H, H1H1 HU nor for a driven system [37| . 
In general, rotational and linear degrees of freedom are correlated in direction [47|. 

In this paper, we investigate the nature of the full energy distribution, that is, the bivariate distribution of rotational 
and translational energy. Motivated by the fact that on average the total energy is not partitioned equally between 
rotational and translational degrees of freedom, we focus on the bivariate distribution P(E, x) of total energy E and 
the modified ratio x = \J E w /E of rotational to total energy. We thereby generalize the understanding of frictionless 
granular matter in terms of the energy distribution to rough grains. 

Our starting point is the nonlinear Boltzmann equation with a collision rule that accounts for the coupling of 
translational and rotational motion due to tangential restitution. We study stationary solutions of the inelastic 
Boltzmann equation that describe steady states achieved through a balance between energy injections that are powerful 
but rare and energy dissipation through inelastic collisions. For high-energy particles we derive a linear equation for 
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the bivariate energy distribution. The latter can be shown to factorize - P(E, x) = p(E)f(x) - into a product of the 
distribution of the total energy, p{E), and the distribution of the fraction of energy stored in the rotational degrees 
of freedom, f(x). The former distribution decays algebraically with energy: p(E) ~ E~ u . The fraction of energy 
stored in rotational motion is universal for energetic particles in the sense that f(x) approaches a limiting distribution 
independent of energy. Furthermore, this quantity has a number of interesting features. First, the distribution is not 
uniform, as it would be, if equipartition were to hold. Second, the distribution is not analytic but has singularities 
at special energy ratios. Third, the distribution and in particular its singularities depend sensitively on the moment 
of inertia and the collision parameters. Only for energetic particles is this distribution well defined. In general, the 
partition of energy into rotational and translational motion depends on the magnitude of the energy. This paper 
specifically addresses two-dimensions, although the theoretical approach and the reported qualitative behavior are 
generic. 

We also develop a general framework for describing high-energy collisions and we use this framework to study freely 
cooling Maxwell Molecules where the moments of the energy distribution can be found in a closed form. For example, 
the two granular temperatures corresponding to the rotational and translational motions are coupled and generally, 
they are not equal. The high-energy behavior found for driven steady-states extends to freely cooling gases. 

The rest of this paper is organized as follows. We review the collision rules and introduce the nonlinear kinetic 
theory in section II. We then derive the linear kinetic theory for high-energy particles in section III. Next, in section 
IV, we study driven steady states and solve for the stationary energy distribution. Freely cooling Maxwell molecules 
are discussed in section V and we conclude in section VI. The Appendices detail technical derivations. 



II. THE NONLINEAR KINETIC THEORY 



Our system consists of an infinite number of identical particles with mass m = 1, radius R, and moment of inertia 
I = qR 2 where < q < 1 is a dimensionless quantity. Each particle has a linear velocity v and an angular velocity 
w. Its total energy is shared by the linear and the rotational motion, E = E v + E w , or explicitly, 

E=\(v* + qRW) (1) 

where v = |v| and w = |w|. 

In a collision between two particles, their velocities (vj , w^) with the labels i = a, 6, change according to 

(v a ,w a ) + (v h ,w h ) -> K,w^) + (vj,wj) (2) 

where the postcollision velocities are denoted by primes. In a binary collision, rotational and translational energy 
are exchanged, while the total energy decreases. In this study, we consider tangential restitution in addition to the 
standard normal restitution. Let be the position of particle i, then the directed unit vector connecting the centers 
of the colliding particles is n = (r — rj,)/\r a — r&|. We term this vector the impact direction. The collision rules are 
most transparent in terms of Uj the particle velocity at the contact point 

u a = v a + i?nxw a (3a) 
u b = v b — R n x w b . (3b) 

The inelastic collision laws state that the normal component of the relative velocity U = u a — is reversed and 
reduced by the normal restitution coefficient < r n < 1. The tangential component is either reversed (rough particles) 
or not (smooth particles) and in any case reduced by the tangential restitution coefficient —1 < r t < 1, according to 
the following collision rules: 

U' • n = -r„U • n, (4a) 
U' x n = -r t U x n. (4b) 

Inelastic collisions conserve linear and angular momentum. Conservation of linear momentum implies that the 
total linear velocity does not change, and conservation of angular momentum enforces that the angular momentum 
of each particle with respect to the point of contact remains the same, because there is no torque acting at the point 
of contact. The collision laws (0} combined with these conservation laws specify the postcollision velocities as linear 
combinations of the precollision velocities [33[ 

fit - Tit 

v = v a — rj n V ■ n n — r\ t (V — V ■ n n) — rj t R n x W w' a =w a -\ n x V H nxnxW (5a) 

qR q 

v 6 = v b + i)nV ■ n n + rjt (V — V ■ n n) + rjtR n x W w' b = wj + — ^ nxV+ — nxnxW (5b) 

qR q 
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where the shorthand notations V = v a — v& and W = w a + wj, were introduced. These collision rules involve the 
normal and tangential collision parameters, defined as 

1 + r„ q 1 + r t 

Vn=—j-, ^d » t = 1 + g 2 ■ (6) 

Their range of values is bounded by 1/2 < 7y„ < 1 and < r/t < q/(l + q). Details of the derivation of the collision 
rules are given in Appendix A, as they are relevant for our discussion. The energy dissipation, AE = E a +E b — E' a — E' b , 
is given by 

A£ = ^(V ■ n) 2 + — ^(V - V ■ nn + iJn x W) 2 . (7) 

4 v ; 1 + q 4 v ' w 

The energy dissipation is always positive, except when the collisions are elastic, r n = 1 and r t = — 1 (perfectly smooth 
spheres) or rt = 1 (perfectly rough spheres). 

The collision rate K(v a , v b ) is the rate by which the two particles approach each other. For hard spheres, this rate 
is simply the normal component of the relative velocity, but we study the general case 

#(v a) v 6 ) = |(v a -v 6 )-n|^ (8) 

with < 7 < 1. Of course, the collision rate vanishes, K — 0, when the particles are moving away from each other, 
(v a — Vb) • n > 0. When particles interact via the central potential r~ K then 7 = 1 — 2^-^ [48]. The two limiting 
cases are hard spheres (7 = 1) and Maxwell molecules (7 = 0) where the collision rate is independent of the velocity 

Giii Hi HI- 

The central quantity in kinetic theory is the probability P(v, w, t) that a particle has the velocities (v, w) at time 
t. We study spatially homogeneous situations where this velocity distribution function is independent of position. 
Under the strong assumption that the velocities of the two colliding particles are completely uncorrelated, the velocity 
distribution obeys the Boltzmann equation 

~~^t W ^ = \ J dil III! dv ad-vr a dv b dw b |(v a - v 6 ) • n| 7 P(v a , w a )P(v 6 ,w 6 ) (9) 
x [<5(v - v^)<5(w - vf' a ) + 5(v ~ v b )<5(w - w b ) - <5(v - v Q )<5(w - w a ) - 6(v - v b )5(w - w&)] . 

We integrate over all impact directions with J dh = 1 and over the precollision velocities weighted by the respective 
probability distributions. There are two gain terms and two loss terms, because the velocities of interest (v, w) can 
be identified with any one of the four velocities in the collision rule and the kernel is simply the collision rate ©. 



III. THE LINEAR KINETIC THEORY 



The focus of this study is the energy distribution that generally depends only on two variables: E v and E w . It 
is our aim to compute the distribution P{E V , E w ) for asymptotically large energies. This will be done for a system 
which is driven at very high energies as well as for an undriven system. 

As a first step to this goal, we simplify the Boltzmann equation in the limit of large energies. Extremely energetic 
particles are rare and as a result it is unlikely that such particles will encounter each other. Hence, energetic particles 
typically collide with much slower particles. Since the collision rules arc linear, the velocity of the slower particle barely 
affects the outcome of the collision. We can therefore neglect the slower velocity. Substituting (v ra , w„) = (vo,wo) 
and (vb,Wb) = (0,0) or (v a ,w a ) = (0,0) and (y b ,w b ) = (vo,Wo) into |5]) gives the cascade process [54 l55l| 

(v ,w ) -> (vi,wi) + (v 2 ,w 2 ) (10) 

where (vq,Wq) is the precollision velocity of the energetic particle and (Vi,Wj) with i = 1,2 are the consequent 
postcollision velocities. With these definitions, the collision rules for extremely energetic particles are 

vi = (1 - ?7„)v • nn+ (1 - %)(v - v • nn) - r? t n X w wi 
V2 = VuVq • fin + ?7 t (vo - v • nn) + rj t h x w w 2 



= U-| 
Vt 

— w 

q 



w 



— n x vo 

q 



Vt - 

— n x v , 

q 



(11a) 
(lib) 
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where we have set R = 1, so that the moment of inertia, I — q, is dimensionless. A collision between a high-energy 
particle and a typical-energy particle produces two energetic particles with an energy total that is smaller than the 
initial energy. This cascade process transfers energy from large scales to small scales. 

Since the cascade process (|TU1) involves only one particle, the tail of the probability distribution P(v, w) obeys the 
linear equation 



<9P(v,w) 
dt 



<intfvo<iwo|vo ■ n| 7 P(vo, wq) [<5(v— vi)5(w— wi) + 5(v — V2)<5(w — W2) — <5(v — vo)<5(w — Wo)l . (12) 



There are two gain terms and one loss term according to the cascade process (|10j) . Formally, this linear rate equation 
can be obtained from the full nonlinear equation © by treating either one of the precollision velocities as negligible 
and then integrating over this small velocity. This procedure leads to four gain terms and two loss terms and thus, 
the factor 1/2 in (J9j) drops out. We stress that the linear equation (I12|) is valid only in the high-energy limit. 

We also comment that the linear equation (|12p for the high-energy tail of the velocity distribution may be valid 
in cases where the full nonlinear equation is not. Whereas the nonlinear equation requires that all possible velocities 
are uncorrelated, the linear equation merely requires that energetic particles are uncorrelated with typical particles. 
This is a much weaker condition. 

In this paper, we restrict ourselves to two space dimensions, i.e. rotating disks. In that case the rotational velocities 
are always perpendicular to the linear velocities. Thus, we conveniently denote the unit vector in the tangential 
direction by t and the unit vector coming out of the plane by z, such that n • t = and n x t = z. The precollision 
velocities of the energetic particle vo = v n h + vtt and wo = wz are compactly written as [v n ,vt,w\. With this 
notation, the postcollision velocities specified in (|12p are 

(1 - T] n )v n , (1 - rj t )v t + r) t w, (r) t /q)v t + (1 - rjt/q)w , and r\ n v n ,r\tivt - w), (r)t/q)(vt - w) , (13) 

respectively. We now treat the three velocity components, namely the normal component of the velocity v n , the 
tangential component of the velocity v t , and the scaled angular velocity y^w as a three dimensional vector with 
magnitude Vo, polar angle 9q, and azimuthal angle 4>o- 

{v n ,Vuy/qw) = ( Vo sin 9 cos 4> , V sin 9 sin 4> , V cos 9 ). (14) 

The magnitude Vq gives the energy Eq = ^V 2 = \{v^ + v 2 + qw 2 ) while the polar angle characterizes the fraction 
of energy stored in the rotational degree of freedom, ^qw 2 /Eq = cos 2 9. In this representation, the postcollision 
velocities are three-dimensional vectors with magnitude Vi, polar angle < 9i < ir, and azimuthal angle < (f>i < 2n. 
The collision rules (fTTj) allow us to express these quantities in terms of V , 9 , 0o : 

(Vi sin^ cos^,T^ sm9i sin<fo, V cos6» 4 ) = (V Ai, VqB,, V Ci) (15) 

where i = 1,2. The magnitudes of the postcollision velocities are proportional to the magnitude of the precollision 
velocity. The three velocity components are scaled by three dimensionless constants Aj, Bi and Cj, that depend on 
the angles 9q> and cf>o of the energetic particle, the collision parameters rj n and rjt, and the moment of inertia q, 





= (1 - rj n ) sin 


9q cos <fio 


(16a) 


Pi 


= (1 - rjt) sin t 


'0 sin (j) + (rft/y/q) cos 9 


(16b) 


Ci 


= (Vt/Vq)sin( 


h sin fa + (1 - rjt/ <?) cos 9o 


(16c) 


A 2 


= r) n sin 9q cos 




(16d) 


P 2 


= Tj t sin 9q sin c, 


^0 - (m/Vo) cosS o 


(16e) 


c 2 


= (Vt/y/q)sm( 


? sin^ - (Vt/q) cos 9 . 


(16f) 



The new energies are proportional to the precollision energies 

E i = a i E , with ai =A 2 + B 2 + C 2 . (17) 

We term the parameters < on < 1 the contraction parameters. Since the collisions are dissipative, these parameters 
satisfy the inequality a\ + a 2 < 1. The equality a% + a 2 = 1 holds only for elastic collisions (r n — \r t \ — 1). The 
energy dissipation is AP = Eq — E\ — E 2 = AP with A = 1 — ot\ — 0,2 or explicitly, 

A = 1 - ? " sin 2 9 cos 2 0o + g 1 ~ - ( sin 2 9 sin 2 O + - cos 2 9 ). (18) 
4 l+a4V a J 
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The polar and azimuthal angles are given by 



cos 0; 



and 



tan < 



Bi 



respectively. 

Let us represent solid angles by = cos 6, (j>. With this definition, the cascade process (fTTj) is 



(19) 



(20) 



with Ei and fi^ given by ([TTJ) and (fT9"|) . Energetic particles have an important property: the solid angle is not coupled 
to the energy! Indeed, the postcollision angles depend only on the precollision angle. The cascade process has the 
following geometric interpretation: a three dimensional vector is duplicated into two vectors. Subsequently, these 
two vectors are scaled down by the contraction parameters (|17l) , and rotated according to the angular transformation 

We can now write the linear Boltzmann equation for P{E, CI), the distribution of energy and solid angle, in a closed 
form 



dP(E,0) 
dt 



= / jdE dn Q yE^smQ cos <j)o\ 1 P(E a ,n )[6(E-E 1 )5(ri-rL 1 ) + 5(E-E2)5(ri-n2)-5(E-E )5(ri-no)]- 



Time was rescaled, t — ► 2 7 / 2 t, to absorb the constant which arises from replacing velocity by energy in the collision 
rate ([8]). Henceforth, we implicitly assume that the distribution P(E, fl) is independent of <j> because the distribution 
of linear velocities must be isotropic. The integration over the energy is performed using the collision rule (|17j) . leading 
to the linear rate equation for the tail of the energy distribution 



dP(E, Sl) 
dt 



— E 1 ! 2 / dO 1 sin 9 cos 4>o I 



P 



E \ <J(fi-fii) 



-P 



Oi2 



1+7/2 



-p(e, n )6{n - n ) 



(21) 



This is a non-local equation as the density of particles with energy E is coupled to the density of particles with the 
higher energies E/a\ and E/a.2- We stress that this equation is a straightforward consequence of the cascade process 
(|20[) and that it can also be derived from the full nonlinear Boltzmann equation. Yet, there may be situations where 
the linear equation (|2"Tj) is valid, while the nonlinear equation ([9]) is not valid. The bivariate energy distributions 
P(E,Cl) and P(E V ,E W ) are completely equivalent but we analyze the former because the cascade process (|2U)) is 
transparent in terms of the total energy and the solid angle. 



IV. DRIVEN STEADY-STATES 



The inelastic Boltzmann equation admits stationary solutions for frictionless particles. These stationary solutions 
describe driven steady-states with rare but powerful injection of energy. The injected energy cascades from high- 
energies down to small energies, thereby balancing the energy lost in collisions. At energies below the injection scale, 
Eqs. ©, (fT2"j) and (f2~Tj) are not altered by the energy source and consequently, the stationary solution of the inelastic 
Boltzmann equation holds up to this large energy scale [13, EH • Here, we seek a corresponding stationary solution 
for particles with rotational degrees of freedom in the high energy limit. 

The stationary solution has to fulfill Eq. (f2"Tj) with the left hand side set to zero 



o sin 6>o cos 0o 



1 



5(n-fii) 



i 



a 



P\- 

1+7/2^ U 2 ' 



ft s(n-n2)-P{E,rio)S{n-n ) 



(22) 



At high-energies, the solid angle is not coupled to the energy, as follows from Eq. (fT9"|) . This fact has a ma- 
jor consequence: the bivariate energy distribution P(E, f2) takes the form of a product of the energy distribution 
p(E) = J dCl P(E,fl) and the distribution of solid angle, g(fi), 



P(E,Sl)^p(E)g(Si) 



(23) 



as E — + oo. The angle distribution is normalized, J dfl g(fl) = 1. It does not depend on the azimuthal angle, because 
on average the two components of the linear velocity are equivalent. Due to the equi-dimensional (in E) structure of 
the steady-state equation ([22]) , the product ansatz (|23|) is a solution when the distribution p(E) decays algebraically 

(24) 



p(E) ~ E-\ 
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FIG. 1: The exponent v for hard spheres (7 = 1) as a function of the coefficients of normal, r n , and tangential, r t , restitution 
coefficients. The numerical procedure for solving (I25|l is detailed below. 



as E — + 00 [5J, [53|. We obtain a closed equation for the distribution g(fi) by substituting the product ansatz |23|) 
with the power-law form into the steady-state equation (j2"2")l 



= / (iilo ff(^o) sin 6*o cos< 



u-l-y/2 



5(f7-^i)+«2 



-1-7/2 



<j(n-n 2 ) -5(n-Q ) 



(25) 



This equation is linear in g(£l). However, it is nonlinear in v and moreover, the solid angles f2,; = fij(Qo) m P9|) and 
the contraction parameters on = a, (i7o) in (| 1 T[) are complicated functions of the solid angle Qq. 

Equation (|25p involves two unknowns quantities, the exponent v and the distribution function g{Q). A solution 
does not exist for arbitrary values of v. In fact, there is one and only one value of v for which there is a solution for 
g(fi). This is the value selected by the cascade dynamics! In other words, (|2l))) is an eigenvalue equation: v is the 
eigenvalue and g(£l) is the eigenfunction. This eigenvalue equation circumvents the full nonlinear equation ([9]) and 
thus, represents a significant simplification. 

The physical interpretation of (|2"5|) involves a cascade process in which the solid angle undergoes a creation- 
annihilation process 



with rate /?o, 
17o — ^ \ Hi with rate 

0,2 with rate /?2- 



(26) 



Here, (3i — | sin0ocos</>o| 7 ai for i = 0,1,2 and ao = 1. There is one annihilation process and two creation processes. 
These processes have relative weights that reflect the powerlaw decay of p(E). At the steady-state, the creation and 
the annihilation terms balance (see Appendix |B|) . as reflected in the integrated form of (|25| 



dflo g(p,o)\ sin#o cos< 



a. 



-1-7/2 



-7/2 



a 



v-l-1/2 



(27) 



< 1. Since o.i < 1, we have the lower bound 



To achieve a steady-state, /3j < /3q for i = 1, 2 and therefore a" x ~ 
v> 1 + 7/2. 

We can immediately check that for elastic collisions, v = 2 + 7/2 (5(| [57| because ai + a.2 = 1, and therefore, we 
conclude the bounds 1 + 7/2 < ^ < 2 + 7/2. The exponent 1^ varies continuously with the restitution coefficients 
r n and r t and the normalized moment of inertia q. This qua ntity must coincide with the value found for frictionless 
particles where tangential restitution is irrelevant (r t — —1) [54. l55l. l58j, but otherwise the exponent is distinct, as 
shown in Fig. [TJ Also, the exponent v increases monotonically with r n and \r t \. We conclude that the rotational 
degrees of motion do affect the power-law behavior (|24|) . 

The azimuthal angle 8 characterizes the fraction of energy stored in the rotational mode, cos 2 9 = E w /E with 
E w = hqw 2 . The angle distribution g(Cl) = (27r) _1 /(cos 9) therefore captures the partition of energy into rotational 
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and translational energies. We introduce the natural variable < x < 1 defined by x = \ cos9\ so that 



x = 




(28) 



and present results for the angle distribution f(x) = 2/(cos#). In equilibrium, energy is partitioned equally into all 
degrees of freedom and therefore <? e q(^) = (47r) _1 or equivalently, 




A. 



Simulation Methods 



We numerically studied the angle distribution f(x) by solving the linear eigenvalue equation (f2"5|) for the "angular" 
process (f26| and by solving the full nonlinear Boltzmann equation ([9]) for the collision process ([2]). Both of these 
equations are solved using Monte Carlo simulations. 

The eigenvalue equation is solved by mimicking the angular process. Throughout the simulation, the value v is 
fixed. There are N particles, each with a given polar angle. A particle with polar angle Qq is picked at random and 
then, a random azimuthal angle 0o is drawn. The polar angles 9\ and 02 are then calculated according to f|19[) . The 
original particle is annihilated with probability (3$ and simultaneously, a new particle with angle 9\ is created with 
probability (3\ and similarly, a second particle with angle 02 is created with probability /?2- Therefore, the number of 
particles may increase by one, remain unchanged, or decrease by one. The exponent v is the value that keeps the total 
number of particles constant in the long time limit. The eigenvalue v is calculated by repeating this simulation for 
various values of v and then using the bisection method [H^] . We present Monte Carlo simulations of 100 independent 
realizations with N = 10 7 particles. 

Driven steady-states are obtained by simulating the two competing processes of inelastic collisions and energy 
injection. In an inelastic collision, two particles are picked at random and also, the impact direction is chosen at 
random. The particle velocities are updated according to the collision law ([5]). Collisions are executed with probability 
proportional to the collision rate. Throughout this process, we keep track of the total energy loss. With a small rate, 
we augment the energy of a randomly selected particle by an amount equal to the loss total and subsequently, reset 
the total energy loss to zero. A fraction of the injected energy is rotational and the complementary fraction is 
translational. We draw this fraction according to the equilibrium distribution (129p . We experimented with different 
angle distributions and found that the resulting stationary state did not change. 

Obtaining the distribution f(x) is generally challenging as it requires excellent statistics. The simulations are 
most efficient for Maxwell molecules because all possible collisions are equally likely. Therefore, for the full nonlinear 
Boltzmann equation @, we present the angle distribution of the energetic particles only for the case 7 = 0. 

For Maxwell molecules, the injection rate is 10 -4 and the system size is N = 10 7 . The corresponding values for 
hard spheres are 10 -2 and N — 10 5 . In all cases, the simulation results represent an average over 10 2 independent 
realizations. Unless noted otherwise, the simulation results are for maximally dissipative (r n — r t — 0) disks (q = 1/2). 



The numerical simulations confirm several of our theoretical predictions. First, the energy distribution approaches a 
steady-state with a power-law high-energy tail. Second, the distribution of the total energy p(E) decays algebraically 
as in (|24p . Third, the exponent v is in excellent agreement with the predictions of the eigenvalue equation. For 
Maxwell molecules, Monte Carlo simulation of the full nonlinear equation yields v = 1.570 ±0.005 whereas numerical 
solution of the eigenvalue equation ([23]) gives v = 1.569 ± 0.005 (Fig. [J). For hard-spheres, where the simulation 
results are slightly less accurate, the corresponding values are v — 2.065 ± 0.005 and v — 2.060 ± 0.005 (Fig. [3]). 
The behavior of the distribution of total energy is therefore qualitatively similar to the behavior in the no-rotation 
case [54l |55|. However, the quantitative behavior is different because the exponent v does depend on the tangential 
restitution coefficient and the moment of inertia (Fig. 



B. The Distribution of Total Energy 



C. The Angle distribution 



The numerical simulations also confirm several of our theoretical predictions concerning the angle distribution. 
Extremely energetic particles have a universal distribution f(x). This distribution is independent of the energy, 
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FIG. 2: The tail of the energy distribution for driven 
Maxwell molecules. Shown are simulation results (solid line) 
and a line with the slope predicted by the theory (dashed 
line). The energy is normalized by the typical energy 1CP 4 . 
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FIG. 4: The angle distribution f(x) obtained by Monte 
Carlo simulation of the angular process (|26[) (solid line) and 
the collision process ((2} (dashed line) for Maxwell molecules. 
The special values xi, X2, and xz discusses in the text are 
indicated by arrows. 
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FIG. 3: The tail of the energy distribution for driven hard 
spheres. Shown are simulation results (solid line) and a line 
with the slope predicted by the theory (dashed line). 
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FIG. 5: The angle distribution f(x) for various collision 
parameters (r n and r t ) for Maxwell molecules. 



provided that the energy is sufficiently large. We had to probe only the most energetic particle out of roughly 10 3 
particles to measure this distribution. For this reason, the linear analysis and the resulting eigenvalue equation are 
valuable because they allow for an accurate and efficient determination of the angle distribution of the energetic 
particles. We also verified that the distribution f(x) obeys the eigenvalue equation (|^5]l . as demonstrated in Fig. 2J 
where the simulations are compared to the solution of the angular process. 

The distribution f(x) has several noteworthy features. First, it is not uniform, implying the breakdown of en- 
ergy equipartition in a granular gas. Furthermore, this distribution is nonanalytic. It contains singularities and 
discontinuous derivatives. There are notable peaks in the distribution so that special values x and special ratios 
E w /E are strongly preferred. The reason for these peaks is the fact that the polar angle is limited. For example, 
cos 2 9 2 < 1/(1 + 0') as seen by substituting cos 6q = ±1 into (fTBjl and (fT§]) . Consequently, there is a special ratio 



Xi 

with the corresponding special energy ratio E w /E 



(30) 



1 + 3 

x\. This is the most pronounced peak in Fig. 2J 
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FIG. 6: The angle distribution f(x) for hard spheres. 



FIG. 7: The angle distribution faii(x) of all particles for 
Maxwell molecules (solid line). Also shown for reference is 
the uniform equilibrium distribution (broken line). 



xi = \/2/3 = 0.81649. Numerically, we observe that the peak becomes more pronounced as the distribution is mea- 
sured at a finer scale, indicating that the distribution function diverges at this point. 

Similarly, there is another special ratio that corresponds to Q\ when cos6*o = ±1, and unlike (|30|) , this location 
depends on the tangential restitution, 



x 2 



1 - Vt/q 



\/ilf/q+(l~7]t/q) 2 



(31) 



Indeed, there is a barely noticeable cusp at X2 — 
"echo" -singularities. For example, using cos#o = x\ and (f>o 



0.942809. Singularities may induce less pronounced 
7r/2 yields the special ratio 



^3 



l + r? t (l-l/g) 



y/q[l- Vt (l -l/qW + {l+Vt(l-l/qW 



(32) 



There is a noticeable peak at the corresponding value 2:3 = -y/50/99 = 0.710669 in Fig. [4j We anticipate that as 
the transformation (|19p is iterated, the strength of the singularities weakens and as a result there are discontinuous 
derivatives of increasing order, a subtle behavior that is difficult to measure. 

The location of the singularities varies with the collision parameters r n and r t and the moment of inertia q. In 
fact, the angle distribution is extremely sensitive to material properties as its shape changes dramatically with these 
parameters, see Fig. [5] The angle distribution also depends on the collision rate and it is much smoother for hard 
spheres, see Fig. [6] Since the collision rate vanishes for grazing collisions, cf> — n/2, the associated singularities 
including in particular (132j) are suppressed. Nevertheless, there is a pronounced jump at the special ratio given by 
([30]) and there are also noticeable cusps. 

The angle distribution of all particles fai\(x) cx J dE P(E,fl) is shown in Fig. [JJ It is substantially different from 
f(x). Therefore, the energy distribution P(E, ft) does not factorize in general and there are correlations between the 
solid angle and the total energy. Only for energetic particles does (l23l) hold. Moreover, f a \i(x) is much smoother in 
comparison with f(x) although there is a jump in the first derivative at the special ratio (|30[) showing that the angle 
distribution of all particles is also non-analytic, see Fig. [7] Generally, the angle distribution depends on energy and 
the deviation from a uniform distribution grows with energy 

We also comment that lone measurement of the moment (a; 2 ) can be misleading. The angle distribution may very 
well have a value close to the equipartition value (x 2 } eq = 1/3 but still, be very far from the equilibrium distribution. 
Indeed, in Fig. [4] (x 2 ) = 0.318, a value that barely differs from the equilibrium value, even though the corresponding 
distribution is far from uniform. The second moment may also differ substantially from the equipartition value and 
for example, (x 2 ) — 0.202 when r n = 0.9 and r t = (Fig. [5]). 

We argue that the qualitative features of the angle distribution should be generic in granular materials. Colli- 
sions involving energetic particles must follow the linear cascade rules (|20l) with the angular transformations (| 19|) . 
The singularities are a direct consequence of these transformations and therefore should be generic. Measuring the 
parameter-sensitive distribution f(x) experimentally is challenging because a huge number of particles must be probed 
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and the measurement has to be accurate. The distribution / a ii(a^) provides a detailed probe of the partition of energy 
into rotational and translation motion. 

V. FREE COOLING 

We now consider freely cooling granular gases that evolve via purely collisional dynamics. Without energy input, 
all energy is eventuall y d issipated and the particles come to rest. This system has been studied extensively [l| for 
hard spheres with (33l.l47| and without rotation [60l |. 

We consider Maxwell molecules where in the absence of rotation an exact treatment is possible [HI IH, H3, EH, E3] • 
When 7 = the Boltzmann equation ^ simplifies 

~~Bt W ^ = \ j dil JJJJ dVad ™ adVbd ™ bP ( Va ' yVa ) p ( Vb >™ b ) ( 33 ) 

x [5(v - v' a )5(w - w' a ) + <S(v - v' b )5(w - w' b ) - S(v - v a )<5(w - w a ) - <5(v - v b )<5(w - w h )]. 
Consequently, the equations for the moments (v n w m ) = JJ dvdxvP(v,w)v n w m close. 

A. The Temperatures 

Here, we consider only the translational temperature defined as the average translational energy, T v — (E v ), and 
the rotational temperature, defined as the average rotational energy T w — (E w ). These two temperatures are coupled 
through the linear equation 



d f T v \ / X vv Xvw \ ( Tv 

dt V T w J \ X wv X W w ) \ T-w 



(34) 



Appendix [C] details the derivation of the matrix of coefficients 

X vv = »7«(1 - Vn) + %(1 ~ Vt) (35a) 
X vw = -2rf/q (35b) 

X wv = ~Vt/ < l (35c) 
X ww = 2(r h /q)(l-r h /q). (35d) 

The two temperatures are coupled as long as i] t ^ or alternatively, r t 7^ —1. 
The solution of is a linear combination of the two eigenvectors 

T w )= C -{L)e- X - t + C + (l)^ t (36) 
with the constants C_ and C+ set by the initial conditions, and c± = (X± — X vv )/X vw . The eigenvalues are 



Xyy ~\~ X W yj / ( A yy A }j 



A±= — f^ ±m — 2 ) +A^A^. (37) 



The larger eigenvalue is irrelevant in the long time limit and therefore, 

T W )- C -{c-) e - Xt ^ 

such that both temperatures decay with the same rate A = A_. Of course, the total temperature also follows the 
same exponential decay, T — T v + T w ~ e~ xt . In this regime, the fraction of rotational energy is on average 

lim ^ = -i=- = - A ~ A "" . (39) 
t^oo T 1 + c_ A + X vw ~ X vv 

The approach toward this value is exponentially fast and the relaxation time is inversely proportional to the difference 
in eigenvalues r = 1/(A + — A_). 

In equilibrium, T w /T =1/3 but for nonequilibrium granular gases the ratio varies. In Fig. 8 we plot the ratio of 
the average rotational energy to the total energy as a function of the coefficients of restitution. In accordance with 
our findings for driven steady-states, energy is not partitioned equally between all the degrees of freedom. 
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FIG. 8: The ratio of average rotational energy to total energy 
as a function of the coefficients of normal, r n , and tangential, 
rt, restitution. 



FIG. 9: The scaling function underlying the energy distri- 
bution (solid line). The distribution was obtained using a 
Monte Carlo simulation with N = 10 7 particles. A dashed 
line with the slope predicted by the theory is also shown for 
reference. 



B. The Energy Distribution 



To study the full energy distribution, it is again convenient to make a transformation of variables from the velocity 
pair (v, w) to the total energy and the solid angle (E, ft). The energy distribution is now time dependent and assuming 
that the temperature - T ~ e~ xt - is the characteristic energy scale we postulate the self-similar form 



P(E,fl,t) -> e xt $(Ee xt ,fl) 



(40) 



with the prefactor ensuring proper normalization, ff dz dft $(2;, ft) — 1. We focus on the high-energy behavior where 
the linear equation (f2Tj) holds. By substituting the scaling form (140)) into this linear equation and setting 7 = 0, we 
find the integro-differential equation governing the scaling function 



\$(z,ft) + \Z — 

dz 



dflo 



1 / z 
— $ — ,fio 



1 / z 

5(fl-fl 1 ) + — $ — ,n 

a 2 \a 2 



s(n - n 2 ) - $(2, n )s(ci - n ) 



.(41) 



We again write the multivariate energy distribution as a product $(z, f2) — > tp(z)g(fl) of the distribution of the total 
energy ip(z) = J dVt $(z, ft) and the distribution of the solid angle g(ft). This form is a solution of the equi-dimensional 
equation (I4ip when the distribution of the total energy decays as a power-law 



at large energies, 



ip(z) - z~ v 

> 00. The angle distribution satisfies the eigenvalue equation 
= y dO . 9 (f7 ){ap 1 (5(fi-O 1 )+a^ 1 (5(fi-^ 2 )- [1 - X(v> - l)]5(fi - fl )}- 



(42) 



(43) 



Of course, setting A = 0, one recovers the steady-state equation (|25|) reflecting that the similarity solution is stationary. 
The factor 1 is replaced by the smaller factor 1 — X(y — 1) that accounts for the constant decrease in the number 
of particles at any given energy because of dissipation. Again, we have a nonlinear eigenvalue equation with the 
eigenvalue v and the eigenfunction gift). 

We solve this eigenvalue equation by performing a Monte Carlo simulation of the same angular process as described 
by but with a different annihilation rate /?o = 1 — \{v — 1). We compare the angle distribution predicted by (|4"3"]) 
with the behavior of the energetic particles in the freely cooling gas. 

The numerical simulations of the inelastic collision process confirm the theoretical predictions. First, the energy 
distribution is self-similar as in (|40|) and the characteristic scale is proportional to the temperature. Second, the 
distribution of the total energy has a power-law tail, as displayed in Fig. and the exponent v is very close to the 
theoretical prediction (numerical simulations of the collision process gives v — 2.98±0.05 while the eigenvalue equation 
yields v = 2.92 ± 0.05). 
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FIG. 10: The angle distribution of the energetic particles. 
Shown are results for the collision process (solid line) and 
for the angular process (dashed line). 




FIG. 11: The angle distribution of all particles for a freely 
cooling gas (solid line). Also shown for reference is the uni- 
form equilibrium distribution. 



The angle distribution deviates even more strongly from the uniform distribution with a very pronounced peak 
(see Fig. ITUf because the dynamics are purely collisional. The singularities are weaker although the one at x\ given 
by (j30|) is clear. The agreement between the solution of the angular process and the Monte Carlo simulations is 
slightly worse than for driven systems because the statistics become prohibitive: now it is necessary to probe the most 
energetic out of roughly 10 6 particles to obtain the asymptotic angle distribution! The sharper power-law decay is 
responsible for this three order of magnitude increase: the cumulative distribution of total energy decays according 
to J E dE'p(E') ~ E~^ with \i = v — 1 about three times larger than before. Finally, the angle distribution of all 
particles deviates only slightly from a uniform distribution (see Fig. ITTj) . We conclude that the behavior of the freely 
cooling gas is qualitatively similar to that found in driven steady-states. 



VI. CONCLUSIONS AND OUTLOOK 



The complete description of granular media with translational and rotational degrees of freedom requires the full 
bivariate distribution of energies. It is not sufficient to consider only the average kinetic energy of translations and 
rotations. Instead the full bivariate distribution is highly nontrivial. We have shown that in the limit of large particle 
energy, this distribution obeys a linear equation. Its solution can be written as a product of two distributions, one 
for the total energy, E = E v + E w , and one for the variable x = \J E w /E, which captures the partition of the 
total energy between rotational and translational motion. The distribution of the total energy decays algebraically 
and the characteristic exponent depends on the collision parameters and the moment of inertia. The variable x is 
not uniformly distributed as in equilibrium. Instead the distribution f(x) is not analytic and displays a series of 
singularities of varying strengths. Remarkably, there are special preferred ratios of rotational-to-total energy. This 
violation of energy equipartition among different degrees of freedom is a direct consequence of the energy dissipation. 
The total energy and the variable x are correlated in general with the deviations from equilibrium increasing with 
energy. These two variable become uncorrelated only at extremely high-energies. 

We have studied both, the system which is driven at extremely high energies and displays a stationary energy 
cascade on energy scales below the driving one, and a freely cooling gas. In the latter gas the bivariate energy 
distribution is time dependent, reflecting the overall decrease of energy. Nevertheless, scaling the total energy with 
temperature, one finds a self-similar form for the distribution, which again factorizes in the high-energy limit. As in 
the driven system, the distribution of the total energy decays as a power law with, however, different exponents for 
the driven and the free cooling system. The angular distribution deviates even more from the uniform (equipartition) 
one in the cooling system. 

It should be straightforward to extend these results to three dimensions where the angular process takes place in 
three dimensions. In the limit of high energies one would again expect a limiting distribution for the partition angle 
x = \J E w jE. Another possible extension refers to a more realistic law of friction, including Coulomb friction (39ll4Q|. 
Finally, it would be of interest to extend the analysis to other systems, where equipartition is violated. An example 
is a binary mixture, where the energy is shared unequally between the two components. 
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APPENDIX A: THE COLLISION RULES 



The total linear momentum ~v' a + ~v' b = v a + Vb is conserved in the collision. The angular momenta of the two 
particles with respect to the point of contact, u:' i7 are given by 

Iu) a = iw a + m R n x v a (Ala) 
IijJb = I~Wb — mRh x Vft. (Alb) 

These are conserved, uj[ = oJi with i = a, b, because there is no torque at the point of contact. In inelastic collisions, 
the normal and tangential components of the relative velocity at the point of contact obey the collision law ((4]) where 
U = V + i?nxW. 

It is convenient to introduce the momentum transfer 8, defined as follows: = v a — 8 and v b = V6+<5. Conservation 
of the angular velocity with respect to the point of contact and Eq. (|A1[) gives — Wi + h x S . In terms of 8, 

the difference in velocity at the point of contact is U' = U — 28 + - n x fi x 8. Substituting this expression into the 
collision laws ([J), the normal and the tangential components of 8 are simply 

8 ■ n = T] n U • n (A2a) 
8 x n = r) t U x n. (A2b) 

Consequently, the momentum transfer is 8 = rj n XJ ■ nn + ^(U — U • nn) or explicitly, 

8 = rin V • fi n + i] t (V - V • n fi) + ij t R n x W (A3) 

We now have the explicit collision rules ((5]). 



APPENDIX B: PARTICLE NUMBER CONSERVATION 



In this appendix, we verify that the stationary solution is consistent with particle number conservation. Maxwell 
Molecules are considered for simplicity. It is straightforward to generalize this calculation to all 7 and to free cooling. 
Our starting point is Eq. (|2"Tj) . specialized to Maxwell molecules, i.e. 7 = 0, 



dP(E,n) 
di 









J dn 











1 

a 2 



-P 



E 

Oil 



As a first step we integrate this equation over the solid angle 

d P (E) 



dt 



CV2 \0!2 



(Bl) 



(B2) 



The power- law behavior (|24[) typically holds in a restricted energy range, £) < E < E u , where £) and E u are upper 
and lower cutoffs. In the driven case, the upper cutoff is set by the energy injection scale. Let N — f^" dEp(E) be 
the total number of particles in this range. With the powerlaw decay ((Ml) , then 

1 



N 



1 



E 



El'") 



(B3) 
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To evaluate this time evolution of N, we substitute the product form ((23)) into (|B2[) and integrate over the energies 
in the aforementioned power-law range, 

^ =Njdn Q g(n ) [a^ + oT 1 -!] • ( B4 ) 
Using Eq. (f27|) . we confirm that the total number of particles is conserved, dN/dt = 0. 

APPENDIX C: THE MATRIX COEFFICIENTS 

In an inelastic collision, the translational energy loss is AE V = E v — E' v with E v = \{v^ + uf) and similarly, the 
rotational energy loss is AE W = E w — E' w with E w = ^q{w 2 + vol). We can conveniently calculate these quantities by 
using v' a = v Q — (5, = Vf, + S, and = w,; + (l/qR)xv x 8, and by expressing the momentum transfer S using the 
natural coordinate system, S — r\ n V n n + r)t(Vt — W)t, 

AE V = r) n (l - Vn )V* + rk(l ~ Vt)V t 2 - V 2 W 2 + %(2 % - l)V t W (Cla) 
AE W = ~(vt/q)V t 2 + 77 t (l - Vt/q)W 2 - Tfe(l - 2r h /q)V t W. (Clb) 

The rate of change of the respective temperatures equals 1/2 the average of this quantities, 4iT v — h{AE v ) and 
-jf.T w = 1;(AE W ). This is seen by multiplying (|33[) by \v 2 and by integrating over the velocity. The averaging is with 
respect to the probability distribution functions of the two colliding particles. The cross-term vanishes, (VtW) = 0, 
by symmetry. Using (V 2 ) = 2{v„) = (v 2 ) = 2T V and (W 2 ) = 2(w 2 ) = AT w /q we obtain the matrix elements (|35]) . 



